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(continued from Part I) 

O 

qth; 3 I/O Monotone systems 

i-H ■ We next describe recent work on monotone input/output systems ("MIOS" from now on). Monotone i/o sys- 
tems originated in the analysis of mitogen- activated protein kinase cascades and other cell signaling networks, 
^ but later proved useful in the study of a broad variety of other biological models. Their surprising breath of 
applicability notwithstanding, of course MIOS constitute a restricted class of models, especially when seen in 
the context of large biochemical networks. Indeed, the original motivation for introducing MIOS, in the 2003 
\Q paper [1], was to study an existing non-monotone model of negative feedback in MAPK cascades. The key 
breakthrough was the realization that this example, and, as it turned out, many others, can be profitably stud- 
ied by decompositions into MIOS. In other words, a non-monotone system is viewed as an interconnection of 
monotone subsystems. Based on the architecture of the interconnections between the subsystems ("network 
structure"), one deduces properties of the original, non-monotone, system. (Later work, starting with [2], 
showed that even monotone systems can be usefully studied through this decomposition-based approach.) 



We review the basic notion from [1]. (For concreteness, we make definitions for systems of ordinary differential 
equations, but similar definitions can be given for abstract dynamical systems, including in particular reaction- 
diffusion partial differential equations and delay-differential systems, see e.g. [3].) The basic setup is that of 
an input/output system in the sense of mathematical systems and control theory [4], that is, sets of equations 

^ =f(x,u), y = h(x), (1) 

in which states x(t) evolve on some subset IC1 W , and input and output values u(t) and y(t) belong to subsets 
U C R m and Y C MP respectively. The coordinates xi, . . . ,x n of states typically represent concentrations 
of chemical species, such as proteins, mRNA, or metabolites. The input variables, which can be seen as 
controls, forcing functions, or external signals, act as stimuli. Output variables can be thought of as describing 
responses, such as movement, or as measurements provided by biological reporter devices like GFP that allow 
a partial read-out of the system state vector (xi, . . . ,x n ). The maps / : X x U — ► W 1 and h : X — > Y are 
taken to be continuously differentiable. (Much less can be assumed for many results, so long as local existence 
and uniqueness of solutions is guaranteed.) An input is a signal u : [0, oo) — > U which is locally essentially 
compact (meaning that images of restrictions to finite intervals are compact), and we write tp(t, xq, u) for the 
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solution of the initial value problem dx/dt(t) = f(x(t),u(t)) with x(0) = xq, or just x{t) if xq and u are clear 
from the context, and y(t) — h(x(t)). See [4] for more on i/o systems. For simplicity of exposition, we make 
the blanket assumption that solutions do not blow-up on finite time, so x(t) (and y(t)) are defined for all 
t > 0. (In biological problems, almost always conservation laws and/or boundedness of vector fields insure 
this property. In any event, extensions to local semiflows are possible as well.) 

Given three partial orders on X, U, Y (we use the same symbol -< for all three orders), a monotone I/O system 
(MIOS), with respect to these partial orders, is a system (pQ) such that h is a monotone map (it preserves 
order) and: for all initial states x±,X2 for all inputs u\,U2, the following property holds: if x\<X2 and u\<U2 
(meaning that u\{t)<U2{t) for all £>0), then ip(t, x\, u)^(p(t, x<i, 112) for all t > 0. Here we consider partial 
orders induced by closed proper cones K C R^, in the sense that x^yiSy — x^K. The cones K are 
assumed to have a nonempty interior and are pointed, i.e. Kf^—K = {0}. A strongly monotone system 
is one which satisfies the following stronger property: if x\ -< X2 and u\ -< U2, then the strict inequality 
(p(t, xi, u) -« (p(t, X2, U2) holds for all t > 0, where x -« y means that y — x is in the interior of the cone K. 

The most interesting particular case is that in which K is an orthant cone in R n , i.e. a set S £ of the form 
{x G l n I SiXi > 0}, where S{ = ±1 for each i. 

When there are no inputs nor outputs, the definition of monotone systems reduces to the classical one of 
monotone dynamical systems studied by Hirsch, Smith, and others [5]. This is what we discussed earlier, 
for the case of orthant cones. When there are no inputs, strongly monotone classical systems have especially 
nice dynamics. Not only is chaotic or other irregular behavior ruled out, but, in fact, almost all bounded 
trajectories converge to the set of steady states (Hirsch's generic convergence theorem [6,7]). 

A useful test for monotonicity with respect to orthant cones, which generalizes Kamke's condition to the i/o 
case, is as follows. Let us assume that all the partial derivatives v>) for i 7^ j, j^(x,u) for all z, j, and 

^(x) for all i,j (subscripts indicate components) do not change sign, i.e., they are either always > or 
always < 0. We also assume that X is convex (much less is needed.) We then associate a directed graph G to 
the given MIOS, with n + m +p nodes, and edges labeled "+" or "— " (or ±1), whose labels are determined 
by the signs of the appropriate partial derivatives (ignoring diagonal elements of df/dx). One may define 
in an obvious manner undirected loops in G, and the parity of a loop is defined by multiplication of signs 
along the loop. (See e.g. [2, 8] for more details.) Then, it is easy to show that a system is monotone with 
respect to some orthant cones in X, U, Y if and only if there are no negative loops in G. A sufficient condition 
for strong monotonicity is that, in addition to monotonicity, the partial Jacobians of / with respect to x 
should be everywhere irreducible . ( "Almost-everywhere" often suffices; see [5,9]. See these references also 
for extensions to non-orthant cones in the case of no inputs and outputs, based on work of Schneider and 
Vidyasagar, Volkmann, and others [10-13]). 

In inhibitory feedback, a chemical species xj typically affects the rate of formation of another species X{ 
through a term like h(xj) = V/(K + xj). The decreasing function h(xj) can be seen as the output of an 
anti-monotone system, i.e. a system which satisfies the conditions for monotonicity, except that the output 
map reverses order: x\ -< X2 => h{x2) di h{x\). 

An interconnection of monotone subsystems, that is to say, an entire system made up of monotone components, 
may or may not be monotone: "positive feedback" (in a sense that can be made precise) preserves monotonicity, 
while "negative feedback" destroys it. Thus, oscillators such as circadian rhythm generators require negative 
feedback loops in order for periodic orbits to arise, and hence are not themselves monotone systems, although 
they can be decomposed into monotone subsystems (cf. [14]). A rich theory is beginning to arise, characterizing 
the behavior of non-monotone interconnections. For example, [1] shows how to preserve convergence to 
equilibria; see also the follow-up papers [3,15-18]. Even for monotone interconnections, the decomposition 
approach is very useful, as it permits locating and characterizing the stability of steady states based upon 
input/output behaviors of components, as described in [2]; see also the follow-up papers [19-21]. 
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Moreover, a key point brought up in [1,2,22,23] is that new techniques for monotone systems in many situations 
allow one to characterize the behavior of an entire system, based upon the "qualitative" knowledge represented 
by general network topology and the inhibitory or activating character of interconnections, combined with 
only a relatively small amount of quantitative data. The latter data may consist of steady-state responses of 
components (dose-response curves and so forth), and there is no need to know the precise form of dynamics or 
parameters such as kinetic constants in order to obtain global stability conclusions and study global bifurcation 
behavior. We now discuss these issues. 

Characteristics 

The main results in [1,2] were built around the study of characteristics, also called step-input steady-state 
responses or (nonlinear) DC gains. To explain this concept, we study the effect of a constant input u(t) = 
uo,t > (in biological terms, the constant input may represent the extracellular concentration of a ligand 
during a particular experiment, for example). For each such constant input, we study the dynamical system 
dx/dt — f(x,uo). Let us assume that all the solutions of this system approach steady states, and let us call 
K(uo) the set of steady states that arises in this way. To each state x in this set K(uq), one may associate the 
corresponding output or measured quantity h(xo). Let k(uo) be the set of all output values that arise in this 
manner. The graph of the set- valued mapping uq i— ► k(uo) is a subset of the cross product space R m x 
which may be though of as a curve when m = p = 1, and which describes the possible steady state output 
values for any given constant input. 

Although many results may be given in more generality, we will assume for the remainder of this paper that 
these mappings are single- valued, not set-valued, in other words that the system is monostable. Thus, a 
(single- valued) characteristic is said to exist for the system if there is a unique steady state for the dynamical 
system dx/dt = f(x,uo), denoted K{u$), and this property is true for all possible constant levels u$. We 
then define the (output) characteristic k : U — > Y as the composition ho K. Under reasonable assumptions 
on X and boundedness, appealing to results from [24, 25] allows one to conclude that K{u$) is in fact a 
globally asymptotically stable ("GAS" from now on) state for dx/dt — f(x,uo), so that all trajectories (for 
this "frozen" value of the input u), converge to K{u$), and the output y(t) converges to k(uo). 

Characteristics (dose response curves, activity plots, steady-state expression of a gene in response to an 
external ligand, etc.) are frequently available from experimental data, especially in molecular biology and 
pharmacology (for instance, in the modeling of receptor-ligand interactions [26]). A goal of MIOS analysis 
is to combine the numerical information provided by characteristics with the qualitative information given by 
(signed) network topology in order to predict global bifurcation behavior. (See [23] for a longer discussion of 
this "qualitative-quantitative approach" to systems biology.) On the other hand, characteristics are also a 
very powerful tool for the purely mathematical analysis of existing models, as we show below. Monotone 
systems with well-defined characteristics constitute a very well-behaved set of building blocks for arbitrary 
systems. In particular, cascades of such systems inherit the same properties (monotone, monostable response). 
The original theorems, in the works [1,2], dealt with systems obtained by interconnecting monotone (or anti- 
monotone) I/O systems with characteristics in feedback. Let us review them next. 

Positive feedback 

The first basic theorem refers to a feedback interconnection of two MIOS 



X\ 

dt 

X2 

dt 



/i(xi,^i), yi = hi(xi) 



(2) 



(3) 



with characteristics denoted by "fc" and "g" respectively. (A degenerate case, in which the second system is 
memory-free, that is, there are no state variables X2 and 1/2 is simply a static function 2/2 (£) — 9( u 2(t)), is also 
allowed. In fact, the proof of the general case can be reduced to that of the degenerate case, simply by taking 
the first system as a cascade connection of the two systems.) 

As in [2], we suppose that the inputs and outputs of both systems are scalar: m\—m2—p\—P2— 1 (see [20] 
for a generalization to high-dimensional inputs and outputs). The "positive feedback interconnection" of the 
systems (Ej) and (j3j) is defined by letting the output of each of them serve as the input of the other (^ 2 =yi="y" 
and u\=y2="u"), as depicted in Figure [D(a). Such positive feedback systems may easily be multi-stable, even 




Figure 1: (a) Positive feedback and (b) characteristics; (c) negative feedback and (d) characteristics 

if the constituent pieces are monostable [27-30]. Let us first discuss how multi-stability may arise in a very 
intuitive and simple example, and later present the general theorem. 

Two typical steady-state responses are as follows. Suppose that P is a protein with Michaelis-Menten pro- 
duction rate and linear degradation: dp/dt = V max u/{k m + u) — kp, where u represents the concentration 
a substrate that is used in P's formation. The reporter variable is y{t)—p{t). In this case, the steady state 
when u(t)=UQ is po—k{uo) — (V rnax /k)uo/(k rn + uq), and we obtain a hyperbolic response, Figl2](a). The re- 




Figure 2: (a) hyperbolic and (b) sigmoidal responses; (c) intersections with degradation; 



sponse in this example is graded ("light-dimmer"): it is proportional to the parameter uo on a large range 
of values until saturation. In contrast, a sigmoidal ("doorbell") response, Figl2](b), which may arise from 
dp/dt = V max u r /(k 1 m + u r ) — kp with Hill coefficient (cooperativity index) r>l, implies that values of uo 
under some threshold will not result in an appreciable activity, while values over the threshold give an abrupt 
change (respectively, po w in steady state and po « V max /k in steady state). Sigmoidal responses are 
critical e.g. if a cell must decide in a binary fashion whether a gene should be transcribed or not, depending 
on an extracellular signal [31-37,37-39,39-45]. Cascades of enzymatic reactions can be made to display such 
ultrasensitive response, if there is a Hill coefficient r>l at each step [46]. 

Multiple attractors may appear if the output y (for example y—p in the example) is fed-back as input u. The 
mechanism might be an autocatalytic process {u—y, e.g. if p helps promote its own transcription) or via a more 
complicated positive feedback pathway from p to u. Formally, substituting u—p into dp/dt = (V rnax u r ) / {k r m + 
u r )—kp (where r—1 or r>l), we obtain the closed-loop equation dp/dt — iy max p r ) / {k r m + p r ) — kp. We plot 
in Fig|2](c) both the first term (formation rate) and the second one (degradation), in cases where r—1 (left) 
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or r>l (right). For r=l, for small p the formation rate is larger than the degradation rate but for large p 
the opposite holds, so the concentration p(t) converges to a unique intermediate value. For r>l, for small p 
the degradation rate is larger, so p(t) converges to a low value, but for large p the formation rate is larger 
and p(t) converges to a high value instead. Thus, two stable states are created, one low and one high, by this 
interaction of formation and degradation. (There is also an intermediate, unstable state.) This reasoning is 
totally elementary, but it provides an intuition for the general result in [2], shown next, which represents a 
far-reaching generalization. (The result can also be viewed as generalizing aspects of the papers [5,47-55], to 
arbitrary MIOS.) 

We consider Fig. H](b), where we have plotted together k and the inverse of g. It is quite obvious that 
there is a bijective correspondence between the steady states of the feedback system and the intersection 
points of the two graphs. With some mild technical conditions of transversality and "controllability" and 
"observability" (the recent papers [56, 57] show that even these mild conditions can be largely dispensed 
with), the following much less obvious facts are true. We first attach labels to the intersection points between 
the two graphs as follows: a label S (respectively, U) if the slope of k is smaller (respectively, larger) than 
the slope of g~ x at the intersection point. One can then conclude that "almost all" (in a measure-theoretic 
sense or in a Baire-category sense) bounded solutions of the feedback system must converge to one of the 
steady states corresponding to intersection points labeled with an S. The proof reduces ultimately to an 
application of Hisrch's generic convergence theorem to the closed-loop system (the technical conditions insure 
strong monotonicity) . However, the value-added is in the fact that stable states can be identified merely from 
the one- dimensional plot shown in Fig. \V[b). (If each subsystem would have dimension just one, one can 
also interpret the result in terms of a simple nullcline analysis; see the Supplementary Section of [19].) We 
remark that the theorems remain true even if arbitrary delays are allowed in the feedback loop and/or if space- 
dependent models are considered and diffusion is allowed (see [23] for a discussion). A new approach [58], 
based not on monotone theory but on a notion of "counterclockwise dynamics," extends in a different direction 
the range of applicability of this methodology. 

We wish to emphasize the potential practical relevance of this result (and others such as [58]). The equations 
describing each of the systems are often poorly, or not at all, known. But, as long as we can assume that 
each subsystem is monotone and uni-stable, we can use the information from the planar plots in Fig. [U(b) 
to completely understand the dynamics of the closed-loop system, no matter how large the number of state 
variables. It is often said that the field of molecular systems biology is characterized by a data-rich/data-poor 
paradox: while on the one hand a huge amount of qualitative network (schematic modeling) knowledge is 
available for signaling, metabolic, and gene regulatory networks, on the other hand little of this knowledge 
is quantitative, at least at the level of precision demanded by most mathematical tools of analysis. On the 
other hand, input/output steady state data (from a signal such as a ligand, to a reporter variable such as 
the expression of a gene monitored by GFP, or the activity of a protein measured by a Western blot) is 
frequently available. The problem of exploiting qualitative knowledge, and effectively integrating relatively 
sparse quantitative data, is among the most challenging issues confronting systems biology. The MIOS 
approach provides one way to combine these two types of data. (For further discussion of this "data-rich/data- 
poor" issue, see [22,23].) 

The theorem from [2] and its generalizations, as well as the negative feedback result discussed below, amount 
to a model-reduction approach. The bifurcation behavior of the complete closed-loop system is obtained from 
a low-order reduction (just to two one-dimensional systems, connected in feedback, when m — p — 1) and 
information on the i/o behavior of the components. This model-reduction view is further developed in [20]. 

More discussion through an example: MAPK cascades 

Mitogen- Activated Protein Kinase (MAPK) cascades are a ubiquitous "signaling module" in eukaryotes, in- 



volved in proliferation, differentiation, development, movement, apoptosis, and other processes [59-61]. There 
are several such cascades, sharing the property of being composed of a cascade of three kinases. The basic 
rule is that two proteins, called generically MAPK and MAPKK (the last K is for "kinase of MAPK," which 
is itself a kinase), are active when doubly phosphorylated, and MAPKK phosphorylates MAPK when ac- 
tive. Similarly, a kinase of MAPKK, MAPKKK, is active when phosphorylated. A phosphatase, which acts 
constitutively (that is, by default it is always active) reverses the phosphorylation. The biological model 
from [19,59] is in FigEfb), were we wrote Zi(t),i = 1,2,3 for MAPK, MAPK-P, and MAPK-PP concentrations 




Figure 3: (a) MAPK cascades; (b) graph; (c) characteristic 



and similarly for the other variables. The input represents an external signal to this subsystem (typically, the 
concentration of a kinase driving forward the reaction). 

We make here the simplest assumptions about the dynamics, amounting basically to a quasi-steady state 
appproximation of enzyme kinetics. (For related results using more realistic, mass-action, models, see [62-64].) 
For example, take the reaction shown in the square in Figl3](a). As (MAPKK-PP) facilitates the conversion 
of z\ into z 2 (MAPK to MAPK-P), the rate of change dz2/dt should include a term a(zi,ys) (and dz\/dt has 
a term —a{z\^y 3 )) for some (otherwise unknown) function a such that a(0,y 3 ) = and > 0, J^- > 
when z\ > 0. (Nothing happens if there is no substrate, but more enzyme or more substrate results in 
a faster reaction.) There will also be a term +/3(z 2 ) to reflect the phosphatase action. Similarly for the 
other species. The system as given would be represented by a set of seven ordinary differential equations (or 
reaction-diffusion PDE's, if spatial localization is of interest, or delay-differential equations, if appropriate). 

This system is not monotone (at least with respect to any orthant cone), as is easy to verify graphically. 
However, as with many other examples of biochemical networks, the system is "monotone in disguise", so to 
speak, in the sense that a judicious change of variables allows one to apply MIOS tools. (Far more subtle forms 
of this argument are key to applications to signaling cascades. A substantial research effort, not reviewed 
here because of lack of space, addresses the search for graph-theoretic conditions that allow one to find such 
"monotone systems in disguise"; see [22,23,63] for references.) 

In this example, which in fact was the one whose study initially led to the definition of MIOS, the following 
conservation laws: yi(i) + y 2 {t) + 2/3 (t) = y tot (total MAPKK) and zi(t) + z 2 {t) + z 3 (t) = z tot (total MAPK) 
hold true, assuming no protein turn-over. This assumption is standard in most of the literature, because tran- 
scription and degradation occur at time scales much slower than signaling. (There is very recent experimental 
data that suggests that turn-over might be fast for some yeast MAPK species. Adding turn-over would lead 
to a different mathematical model.) These conservation laws allow us to eliminate variables. The right trick 
is to eliminate y 2 and z 2 . Once we do this, and write y 2 — y to t ~ V\ ~ V3 an d z 2 = z tot — z\ — 23, we are left 
with the variables x, yi, 2/3, z\, z 3 . For instance, the equations for Zi,z 3 look like: 

^jT = ~ a ( z i,ys) +/?0tot -z 1 -zs) ^ = 7(z tot - zi - z 3l y 3 ) - 8{z 3 ) 

ft 



for appropriate increasing functions a:, /?, 7, 5. The equations for the remaining variables are similar. The 
graph, ignoring, as usual, self-loops (diagonal of Jacobian), is shown in FigJ3](b). This graph has no negative 
undirected loops, showing that the (reduced) system is monotone. A consistent spin assignment (including 
the top input node and the bottom output node) is shown in Figure [U It is also true that this system has a 




Figure 4: Consistent assignment for simple MAPK cascade model 



well-defined monostable state space response (characteristic); there is no space to discuss the proof here, so 
we refer the reader to the original papers [1,8]. 

Positive and negative feedback loops around MAPK cascades have been a topic of interest in the biological 
literature. For example, see [37,40] for positive feedback and [65,66] for negative feedback. Since we know 
that the system is monotone and has a characteristic, MIOS theory as described here can indeed be applied 
to the example. We study next the effect of a positive feedback u = gy obtained by "feeding back" into the 
input a scalar multiple g of the output. (This is a somewhat unrealistic model of feedback, since feedbacks 
act for example by enhancing the activity of a kinase. We pick it merely for illustration of the techniques.) 

The theorem does not require actual equations for its applicability. All that is needed is the knowledge that 
we have a MIOS, and a plot of its characteristic (which, in practice, would be obtained from interpolated 
experimental data). In order to illustrate the conclusions, on the other hand, it is worth discussing a particular 
set of equations. We take equations and parameters from [19,22,23]: 

dx V2 x 

— — ; VVqU + Vx 

at K2 + x 

dyi _ ve (ytot - ?/i - 2/3) _ V3%yi 
dt k 6 + (y tot -yi-ys) k 3 + y 1 

dm = V4x(y tot - yi - y 3 ) _ v$ y 3 
dt fc 4 + (y t ot - yi-ys) h + y 3 

dz\ _ v 10 (z tot -zi-zs) v 7 y 3 zi 
dt fcio + (z tot -Z1-Z3) k 7 + zi 

dz 3 _ v 8 ys (gtot -Z1-Z3) _ v 9 z 3 
dt k 8 + (z tot - Z1-Z3) k 9 + z 3 

with output Z3. Specifically, we will use the following parameters: ^0 = 0.0015, v\ = 0.09, V2 = 1.2, ^3 = 0.064, 
V4 = 0.064, v 5 = 5, v 6 = 5, v 7 = 0.06, v 8 = 0.06, v 9 = 5, v 10 = 5, y tot = 1200, z tot = 300, k 2 = 200, k 3 = 1200, 
k A = 1200, h = 1200, k 6 = 1200, k 7 = 300, k 8 = 300, k 9 = 300, k 10 = 300. (The units are: totals in nM 
(mol/cm 3 ), v's in nM-sec -1 and sec -1 , and fc's in nM.) 

With these choices, the steady state step response is the sigmoidal curve shown in Figj3](c), where y is the 
output £3. We plotted in the same figure the inverse g~ x of the characteristic of the feedback system, in this 
case just the linear mapping y = (l/g)u, for three typical "feedback gains" (#=1/0.98, 1/2.1, 1/6). 
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For g = 1/0.98 (line of slope 0.98 when plotting y against u), there should be a unique stable state, with a high 
value of the output y = and trajectories should generically converge to it. Similarly, for g = 1/2.1 (line of 
slope 2.1) there should be two stable states, one with high and one with low y = z$, with trajectories generically 
converging to one of these two, because the line intersects at three points, corresponding to two stable and one 
unstable state (exactly as in the discussion concerning the simple protein formation/degradation sigmoidal 
example in Figj2](c)). Finally, for g = 1/6 (line of slope 6), only the low-?/ stable state should persist. Figj5](a- 
c) shows plots of the hidden variable ys(t) (MAPKK-PP) for several initial states, confirming the predictions. 
The same convergence results are predicted if there are delays in the feedback loop, or if concentrations depend 




Figure 5(a),(b),(c) y 3 , g = 1/0.98, 1/2.1, 1/6 

on location in a convex spatial domain. Results for reaction-diffusion PDE's and delay-differential systems 
are discussed in [23], and simulation results for this example are also provided there. 

We may plot the steady state value of y, under the feedback u — gy, as the gain g is varied, Fig|6](a). 
This resulting complete bifurcation diagram showing points of saddle-node bifurcation can be also completely 




Figure 6: (a) bifurcation diagram and relaxation (b) oscillation (ys) 

determined just from the characteristic, with no need to know the equations of the system. Relaxation 
oscillations may be expected under such circumstances if a second, slower, feedback loop is used to negatively 
adapt the gain as a function of the output. Reasons of space preclude describing a very general theorem, 
which shows that indeed, relaxation oscillations can be guaranteed in this fashion: see [18] for technical details, 
and [23] for a more informal discussion. Figj6](b) shows a simulation confirming the theoretical prediction 
(details in [18,23]). 

Negative feedback 

A different set of results apply to inhibitory or negative feedback interconnections of two MIOS systems ([2j)- 
([3]). A convenient mathematical way to define "negative feedback" in the context of monotone systems is to 
say that the orders on inputs and outputs are inverted (example: an inhibition term of the form j^_ y , as usual 
in biochemistry). Equivalently, we may incorporate the inhibition into the output of the second system 0, 
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which is then seen as an anti-monotone I/O system, and this is how we proceed from now on. See Fig(T](c). 
We emphasize that the closed-loop systems that result are not monotone, at least with respect to any known 
order. 

The original theorem, from [1], is as follows. We assume once more that inputs and outputs are scalar 
(ra=p=l; see [3] for generalizations). We once again plot together k and g -1 , as shown in FigJTfd). Consider 
the following discrete iteration: 

^+1 = (g°k)(ui). 

Then, if solutions of the closed-loop system are bounded and if this iteration has a globally attractive fixed 
point u, as shown in FigQ](d), then the feedback system has a globally attracting steady state. (An equivalent 
condition, see [3], is that the iteration have no nontrivial period-two orbits.) We call this result a small gain 
theorem ("SGT"), because of its analogy to concepts in control theory. 

It is easy to see that arbitrary delays may be allowed in the feedback loop. In other words, the feedback could 
be of the form u(t) — y{t — h), and such delays (even with h = h(t) time varying or even state-dependent, 
as long as t — h{t) — > oo as t — > oo) do not destroy global stability of the closed loop. In [17], we have now 
shown also that diffusion does not destroy global stability either. In other words, a reaction-diffusion system 
(Neumann boundary conditions) whose reaction can be modeled in the shown feedback configuration, has the 
property that all solutions converge to a (unique) uniform in space solution. This is not immediately obvious, 
since standard parabolic comparison theorems do not immediately apply to the feedback system, which is not 
monotone. 

Example: MAPK cascade with negative feedback. 

As with the positive feedback theorem, an important feature is applicability to highly uncertain systems. As 
long as the component systems are known to be MIOS, the knowledge of I/O response curves and a planar 
analysis are sufficient to conclude GAS of the entire system, which may have an arbitrarily high dimension. 
For example, suppose we take a feedback like u^a+b/ (c+zz), with a graph as shown in FigH(a) 3 which also 
shows the characteristic and a convergent discrete 1-d iteration [23]. Then, we are guaranteed that all solutions 
of the closed-loop system converge to a unique steady state, as confirmed by the simulations in Fig0(b), which 
shows the concentrations of the active forms of the kinases. 
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Figure 7: inhibition: (a) spiderweb and (b) simulation 
Example: testosterone model. 

This example is intended to show that even for a classical mathematical biology model, a very simple appli- 
cation of the result in [1] gives an interesting conclusion. The concentration of testosterone in the blood of 
a healthy human male is known to oscillate periodically with a period of a few hours, in response to similar 
oscillations in the concentrations of the luteinising hormone (LH) secreted by the pituitary gland, and the 
luteinising hormone releasing hormone (LHRH), normally secreted by the hypothalamus (see [67]). The well- 
known textbook [68] (and its previous editions) presents this process as an example of a biological oscillator, 

Q 




and proposes a model to describe it, introducing delays in order to obtain oscillations. (Since the textbook 
was written, the physiological mechanism has been much further elucidated, and this simple model is now 
known not to be correct. However, we want merely to illustrate a point about mathematical analysis.) The 
equations are: 

A 

R = -7? — ™ - hR 
K + T 

L giR — b 2 L 

t = g 2 L(t-r)-b 3 T 

(i?, L,T = concentrations of hormones luteinising hormone releasing, luteinising, and testosterone, r = delay; 
we use "x" to denote time derivative). The system may be seen as the feedback connection of the MIOS 
system 

R = u-biR 
L giR — b 2 L 
f = g 2 L-b 3 T 

with the inhibitory feedback u(t) = g(T — r) = A/(K + T(t — r)) after moving the delay to the loop (without 
loss of generality). The characteristic is linear, T — k{u) — b g J^l 3 so g o k is a fractional transformation 
S( u ) — Since such a transformation has no period- two cycles, global stability follows. (For arbitrary, 

even time-varying, delays.) This contradicts the existence of oscillations claimed in [68] for large enough 
delays. (See [69], which also explains the error in [68].) 

Example: Lac operon. The study of E. Coli lactose metabolism has been a topic of research in mathemat- 
ical biology since Jacob and Monod's classical work which led to their 1995 Nobel Prize. For this example, we 
look at the subsystem modeled in [70]. The lac operon induces production of permease and /3-gal, permease 
makes the cell membrane more permeable to lactose, and genes are activated if lactose present; lactose is 
digested by the enzyme /3-gal, and the other species are degraded at fixed rates. (In this model from [70], 
lactose and isolactose are identified, and catabolic repression by glucose via cAMP is ignored.) Delays arise 
from translation of permease and /3-gal. The equations are: 

x\{t) — g{x±(t — r)) — b\Xi(t) lac operon mRN A 

x 2 {t) — x\(t) — b 2 x 2 {t) /3-galactoside permease 

xs(t) = rx\(t) — 63X3 (t) /3-galactosidase 

x±(t) — Sx 2 (t) — xs(t)x4 (t) lactose 

with g(x) := (1 + Kx p )/{1 + x p ), K > 1, and the Hill exponent p representing a cooperativity effect. (All 
delays have been lumped into one.) We view this system as a negative feedback loop, where u=x\, v—x^ of a 
MIOS system (details in [3]). Since there are two inputs and outputs, now we must study the two-dimensional 
iteration 

g(v) Sbxb^u 
b ' rb 2 g(v) 

Based on results on rational difference equations from [71], one concludes that there are no nontrivial 2- 
periodic orbits, provided that p < {y/~K+l)/{K— 1), for arbitrary 61, b 2 , 63, r, S. Hence, by the theorem, there 
is a unique steady state of the original system, which is GAS, even when arbitrary delays are present, 

These and other conditions are analyzed in [3] , where it is also shown that the results from [70] are recovered 
as a special case. Among other advantages of this approach, besides generalizing the result and giving a 
conceptually simple proof, we have (because of [17]) the additional conclusion that also for the corresponding 
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(u,v) ^ (gok)(u,v) 



reaction-diffusion system, in which localization is taken account of, the same globally stable behavior can be 
guaranteed. 

Example: Circadian oscillator. As a final example of the negative feedback theorem, we pick Goldbeter's 
[72, 73] original model of the molecular mechanism underlying circadian rhythms in Drosophila. (In this 
oversimplified model, only per protein is considered; other players such as tim are ignored.) PER protein 
is synthesized at a rate proportional to its mRNA concentration. Two phosphorylation sites are available, 
and constitutive phosphorylation and dephosphorylation occur with saturation dynamics, at maximum rate 
ViS and with Michaelis constants K{. Doubly phosphorylated PER is degraded, also satisfying saturation 
dynamics (with parameters v^kd), and it is translocated to the nucleus with rate constant k\. Nuclear 
PER inhibits transcription of the per gene, with a Hill-type reaction of cooperativity degree n and threshold 
constant Kj, and mRNA is produced, and translocated to the cytoplasm, at a rate determined by a constant 
v s . Additionally, there is saturated degradation of mRNA (constants v m and k m ). The model is {Pi — per 
phosphorylated at i sites, Pjy = nuclear per, M = per mRNA): 

M = Vs K I /{K I + P^)-v m M/{k m + M) 

Po = k s M-V 1 P /{K 1 +P ) + V 2 P 1 /{K 2 + P 1 ) 

A = V 1 P /{K 1 + P )-V 2 P 1 /{K 2 +P 1 )-V S P 1 /{K S + P 1 ) + V 4 P 2 /{K 4 + P 2 ) 
P 2 = v 3 p 1 /{K 3 + P 1 )-V 4 P 2 /{K 4 + P 2 )-k 1 P 2 + k 2 P N -v d P 2 /{k d + P 2 ) 
Pn = hP 2 -k 2 P N . 

Parameters are chosen exactly as in Goldbeter's original paper, except that the rate v s of mRNA translocation 
to the cytoplasm is taken as a bifurcation parameter. The value v s = 0.76 from [72] gives oscillatory behavior. 
On the other hand, we may break up the system into the M and Pi,Pn subsystems. Each of these can 
be shown to be MIOS and have a characteristic. (The existence of a characteristic for the P-subsystem is 
nontrivial, and involves the application of Smillie's Theorem [74] for strongly monotone tridiagonal systems, 
and more precisely, repeated application of a proof technique in [74] involving "eventually monotonicity" of 
state variables.) When t; 5 =0.4, the discrete iteration is graphically seen to be convergent (see FigJH](a)), so 
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Figure 8: (a) convergent iteration; (b) divergent iteration; (c) oscillations 

the theorem guarantees global asymptotic stability even when arbitrary delays are introduced in the feedback. 
Bifurcation analysis on delay length and v s indicates that local stability will fail for somewhat larger values. 
Using again the graphical test, we observe that for v s =0.5 there appears limit cycle for the discrete iteration 
on characteristics, see FigJH^b). This suggests that oscillations may exist in the full nonlinear differential 
equation, at least for appropriate delays lengths. Indeed, the simulation in Figj8](c) displays such oscillations 
(see [8,14]). 
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Multivalued Characteristics 



For simplicity, we have not discussed the case when characteristics are set- valued instead of single- valued. 
This general case can also be productively studied with the toolkit afforded by MIOS interconnection theory, 
see [56,57] for positive feedback and [21] for negative feedback. 

4 Conclusions 

There is a clear need in systems biology to study robust structures and to develop robust analysis tools. The 
theory of monotone systems provides one such tool. Interesting and nontrivial conclusions can be drawn from 
(signed) network structure alone, which is associated to purely stoichiometric information about the system, 
and ignores fluxes. 

Associating a graph to a given system, we may define spin assignments and consistency, a notion that may 
be interpreted also as non-frustration of Ising spin-glass models. Every species in a monotone system (one 
whose graph is consistent) responds with a consistent sign to perturbations at every other species. This 
property would appear to be desirable in biological networks, and, indeed, there is some evidence suggesting 
the near-monotonicity of some natural networks. Moreover, "near" -monotone systems might be "practically" 
monotone, in the sense of being monotone under disjoint environmental conditions. 

Dynamical behavior of monotone systems is ordered and "non chaotic". Systems close to monotone may be 
decomposed into a small number of monotone subsystems, and such decompositions may be usefully employed 
to study non-monotone dynamics as well as to help detect bifurcations even in monotone systems, based only 
upon sparse numerical data, resulting in a sometimes useful model-reduction approach. 
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